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Abstract. It is well-established that including spatial structure and stochastic noise in 
models for predator-prey interactions invalidates the classical deterministic Lotka-Volterra 
picture of neutral population cycles. In contrast, stochastic models yield long-lived, but 
, ultimately decaying erratic population oscillations, which can be understood through a resonant 

■ amplification mechanism for density fluctuations. In Monte Carlo simulations of spatial 

stochastic predator-prey systems, one observes striking complex spatio-temporal structures. 
These spreading activity fronts induce persistent correlations between predators and prey. In 
• the presence of local particle density restrictions (finite prey carrying capacity), there exists an 

extinction threshold for the predator population. The accompanying continuous non-equilibrium 
phase transition is governed by the directed-percolation universality class. We employ field- 
theoretic methods based on the Doi-Peliti representation of the master equation for stochastic 
particle interaction models to (i) map the ensuing action in the vicinity of the absorbing state 
^ . phase transition to Reggeon field theory, and (ii) to quantitatively address fluctuation-induced 

(N ■ renormalizations of the population oscillation frequency, damping, and diffusion coefficients in 

\l , the species coexistence phase. 
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1. Introduction 

Over the past decade, mathematical and computational tools from statistical physics have 
been increasingly and quite successfully applied to ecological problems, including attempts at a 
^ I quantitative understanding of biodiversity In this context, physicists typically consider 

^ ■ simplified idealized models that hopefully capture the essential features of interacting biosystems; 

- - - leaving aside some of the biological complexity allows the consistent incorporation of stochastic 

fluctuations and spatio-temporal correlations, whose crucial importance has long been recognized 
in the field but is still often neglected. 

Predator-prey models defined via reaction-diffusion systems on a regular lattice, whose 
rate equations in the well-mixed mean-field limit reduce to the classic coupled Lotka-Volterra 
ordinary differential equations, constitute paradigmatic examples of the dynamics of two 
competing populations [6]-[8]. Monte Carlo simulations of these models, specifically in two 
dimensions, display a remarkable wealth of intriguing features (for a fairly recent overview, 
see, e.g., Ref. [9]): In contrast to the regular non-linear oscillations of the deterministic Lotka- 
Volterra model for which the population densities invariably return to their initial values (c.f. 
figure [1] below), computer simulations display persistent, but eventually decaying stochastic 
population oscillations (figure [2]) [IO]-[l7]. In the absence of spatial degrees of freedom, these 
erratic population oscillations may be understood through a resonant stochastic amplification 



mechanism [18] that drasticahy extends the transient time interval before any finite system 
ultimately reaches its absorbing stationary state, where the predator population becomes 
extinct fT9]. In spatially extended systems, it is well-known that the mean-field Lotka-Volterra 
reaction-diffusion equations allow for traveling wave solutions [20]-[22]. In the corresponding 
stochastic spatial realizations, spreading activity fronts (figure [U [23]) induce short-ranged but 
significant positive correlations of either species, and anti-correlations between the predator 
and prey populations, which have the effect of further enhancing the amplitude and life time 
of local population oscillations [23] • We have investigated various different variants of 
stochastic spatial Lotka-Volterra models for competing predator-prey populations, and found 
these intriguing spatio-temporal structures to be remarkably robust against rather drastic 
changes of the detailed microscopic interaction rules [Ml ES]) even the introduction of 
quenched spatial disorder in the reaction rates [2B] . 

In this brief communication, I will provide an overview of our Monte Carlo simulation results, 
specifically contrasting model variants with and without restrictions on the number of particles 
per lattice site. The former describe ecological systems with finite local carrying capacity, and 
display a continuous non-equilibrium phase transition from an active species coexistence state 
to an absorbing phase wherein the predators become extinct. Numerical evidence supports 
the general expectation [27]-[32] that this extinction transition should be governed by the 
directed-percolation universality class [3 [8], [II]-[I1], [111 [T7]. I will then demonstrate how 
field-theoretic tools based on the Doi-Peliti representation of the master equation for stochastic 
interacting particle systems [33]-[35] (for recent reviews, see Refs. [36ll37J), augmented with 
a means to incorporate restricted site occupation numbers [38], can be employed to gain a 
comprehensive understanding of fluctuation and correlation effects in Lotka-Volterra predator 
models. Specifically, the effective action near the extinction transition in model variants with 
restricted site occupations will be explicitly mapped onto Reggeon field theory which describes 
the universal scaling of directed-percolation clusters [271 [32 . 39, 40j. Moreover, expanding on the 
treatment in Ref. [41|, I shall report a computation of the fluctuation-induced renormalizations of 
the population oscillation frequency, damping, and diffusion coefficients in the species coexistence 
phase to lowest order in a perturbation expansion with respect to the predation rate [32] . 

2. Stochastic lattice Lotka Volterra models 
2.1. Model variants and mean-field description 

We consider a two-species system of diffusing particles (with diffusion constant D) that undergo 
the following stochastic reactions: 

^ — with rate 

A + B ^ A + A with rate A, (1) 
B ^ B + B with rate a. 

The 'predators' A decay or die spontaneously at rate /i > 0, whereas the 'prey' B produce 
offspring with rate cr > 0. In the absence of the binary 'predation' interaction with rate A, the 
uncoupled first-order processes would naturally lead to predator extinction a{t) = a(0) e"'^*, and 
Malthusian prey population explosion h{t) = 6(0) e""*; here a{t) and h{t) respectively indicate the 
A / B concentrations or population densities. The binary predation reaction induces species 
coexistence through the non-linear interaction of both particle species. 

In the simplest spatial realization of this stochastic reaction-diffusion model, both particle 
species are represented by unbiased random walkers on a d-dimensional hypercubic lattice, and 
we allow an arbitrary number of particles per lattice site |24j. All reactions ([T|) can then be 
implemented strictly on-site: Offspring particles are placed on the same lattice point as their 
parents, and the predation reaction happens only if an A and a B particle meet on the same 



lattice site. If we then assume the populations to remain well mixed, and consequently ignore 
both spatial fluctuations and correlations, we can approximately describe the coupled reactions 
([1]) through the associated mean-field rate equations for spatially homogeneous concentrations 
a{t) = {a{x,t)), b{t) = {b{x,t)), where a{x,t) and b{x,t) respectively denote the local predator 
and prey densities. One then arrives at the classic Lotka-Volterra equations [3], a coupled set 
of two ordinary non-linear differential equations: 

a{t) = Xa{t)b{t) - fj,a{t) , b{t) = a b{t) - X a{t) b{t) . (2) 



The rate equations ^ display three stationary states {as,bs), namely the empty absorbing 
state with total population extinction (0,0), which is obviously linearly unstable if a > 0; 
a predator extinction absorbing state wherein the prey population diverges (0, oo), which for 
A > is also linearly unstable; and finally a species coexistence state (a„ = a/X,bu = fJ-fX), 
which however represents only a marginally stable fixed point with purely imaginary eigenvalues 
±iy/]Ta of the associated Jacobian stability matrix: Linearizing eqs. ([2]) near {au,bu) results 
in the coupled differential equations 6a{t) = a6b{t), Sb{t) = —ij.6a{t), which are readily solved 
by 6a{t) = 6a{0) cos {^fj, a t) +5b{0) \J a j [i sin (^fiat) and 6b{t) = —6a{0) \J [xja sin (y/JTat) + 
(56(0) cos (^fiat), describing harmonic oscillations about the center fixed point {au,bu) with 
frequency uj = 2tt f = s/JTa. Indeed, the phase space trajectories for the full non- linear coupled 
differential equations ([2]) are determined by da/db = a{Xb — fi)/b{a — Xa), with a conserved 
first integral 

K{t) = X[a{t) + b{t)] -a In a{t) - /i In b{t) = K{0) . (3) 

Consequently, as depicted in figured! the solutions of the deterministic mean-field Lotka-Volterra 
model are closed orbits in phase space, i.e., regular periodic non-linear population oscillations 
whose amplitudes are fixed by the initial configuration. Naturally, the precise periodic return to 
the initial concentration values does not appear to be a very realistic feature. In addition, the 
neutral cycles of the coupled mean-field rate equation system ([2]) indicates that this deterministic 
mathematical model is fundamentally unstable with respect to slight modifications [4J. 
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Figure 1. Solutions of the coupled Lotka-Volterra mean-field rate equations ((2|): (a) Non-linear predator 
a{t) (red) and prey b{t) (blue) population density oscillations; (b) periodic orbits in the a-b phase plane. For 
small amplitudes the oscillations become harmonic (circular orbits) with frequency lo = ^fJTa. (Reproduced with 
permission from Ref. 2 1 . p. 4.) 



One such modification that aims at rendering the Lotka-Volterra system more relevant 
biologicahy is to introduce a finite carrying capacity (total particle density) p > that limits 
the prey population growth, modeling, e.g., the effect of limited food resources [4J. Within the 
mean-field rate equation approximation, the second differential equation in ([2]) is then replaced 
with 

b{t) = a h{t) [1 - b{t) Ip\ - A a{t) b{t) . (4) 

The non-trivial stationary states in this restricted Lotka-Volterra model are predator extinction 
and prey saturation (0, p), linearly stable for X < Xc = p/p; and species coexistence {ar,br) with 
br = p/X and = {1 — p/ pX)a/ X, which both exists and is linearly stable provided the predation 
rate is sufficiently large, A > Ac- The eigenvalues of the Jacobian now acquire negative real parts, 
e± = —pa 1 ± ^/T^^'ApXJpX/]I^^TyJo■ /2pX, which implies an exponential approach to the 
stable fixed point {ar,br), replacing the neutral cycles of the unrestricted model i^. Moreover, 
for a > as = 4:Xp{pX/p — 1) > 0, or p/p < A < Ag = ^1 + ^/l + a/p^ /^/2/o, the eigenvalues 
are real, indicating a nodal stable fixed point, whereas for a < as or X > Xs, i.e., deep in the 
species coexistence phase, the eigenvalues e± turn into a complex conjugate pair, and {ar,br) 
becomes a stable spiral singularity which is approached in a damped oscillatory manner. Adding 
spatial degrees of freedom, finite local carrying capacities can be implemented in a lattice model 
through limiting the maximum occupation number per site for each species. Most drastically, 
one can permit at most a single particle per lattice site [9]; the binary predation reaction then 
has to occur between predators and prey on adjacent nearest-neighbor sites, and new offspring 
needs to be placed on neighboring positions. In that case, one can in fact entirely dispense with 
hopping processes, since all particle production reactions entail population spreading as well. 

In summary, already within the mean-field rate approximation, a finite prey carrying capacity 
p, which can be viewed as the average result of local restrictions on the prey density originating 
from limited resources, crucially changes the phase diagram: There emerges an extinction 
threshold (at Ac for fixed p) for the predator population, which in a spatially extended 
system becomes a genuine continuous active-to-absorbing non-equilibrium phase transition in 
the thermodynamic limit of infinite system size and time. 



2.2. Monte Carlo simulation results 

Various authors have studied stochastic lattice predator-prey models that in the well-mixed 
mean-field limit reduce to the classical Lotka-Volterra system [6]-[8], |10j-|17j. In this section, I 
shall briefiy discuss the pertinent results from our own individual-based Monte Carlo simulation 
studies, performed mostly on two-dimensional square lattices with periodic boundary conditions. 
Technical details and more precise descriptions of the algorithms we have employed can be found 
in Refs. i9j and [24J. 

Figure [2] shows typical simulation data for the temporal evolution of the total predator and 
prey particle densities in a two-dimensional stochastic lattice model with (almost) arbitrarily 
large site occupation numbers and on-site reactions [23]. One observes long-lived but clearly 
damped population oscillations that are actually quite independent of the initial state; neither 
are they caused by the constancy of the first integral K, eq. ([3]), that follows from the 
deterministic rate equations: As is apparent from the numerical data, in the stochastic spatial 
model K{t) is manifestly time-dependent, and in fact traces the overall population oscillations. 
We note that as the system size increases, the relative oscillation amplitudes become smaller; in 
the thermodynamic limit, the quasi-periodic population fluctuations eventually die out entirely. 
From the marked peaks in the Fourier-transformed concentration signals, a{f) = J a{t) e^'^*-^* dt 
for the predators, and similarly for the prey density, we may infer a characteristic oscillation 
frequency /. As illustrated in figure [3l the typical population oscillation frequencies thus 
obtained roughly follow the square-root dependence on the rates p and a as predicted by 
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Figure 2. Monte Carlo simulation data for a stochastic spatial Lotka-Volterra system on a 1024 x 1024 square 
lattice with periodic boundary conditions, in the absence of site occupation number restrictions: (a) Temporal 
evolution for the predator a{t) (red) and prey b{t) (blue) population densities, and the quantity K{t) (green) for 
a — 0.1, n — 0.2,, and A = 1.0; (b) Fourier-transformed population density signals |a(/)| and |6(/)| for a = 0.03, 
H = 0.1,, and A — 1.0. (Reproduced with permission from Ref. [24], pp. 9 and 10.) 

the linearized mean-field approximation, but with measurable deviations both for low and 
high rates. Yet the numerical frequency values are reduced by about a factor of four in 
the stochastic spatially extended system, an apparent considerable downward renormalization 
caused by fluctuations and reaction-induced spatio-temporal correlations [21] . Note also that 
figure [3|^ a) shows a remarkably similar functional dependence of / on the rates /x and a. 




Figure 3. Characteristic peak frequencies obtained form |a(/)| or \b(f)\, (a) as functions of the rates a (red 
squares) and fi (blue diamonds), with otherwise a = 0.1 — fi and A — 1.0 held fixed on a 1024 lattice without site 
occupation restrictions; (b) as function of a with fi = 0.1, A — 1.6, D = on a 128 x 128 square lattice with at 
most a single particle per lattice site (black line) . The dashed black lines represent the oscillation frequency from 
linearized mean-field theory, / — y/JTa. (Reproduced with permission from Ref. [21], p. 10 and Ref. [S], p. 469.) 



Very similar features are found in stochastic spatial Lotka-Volterra models that incorporate 
stringent site occupation number restrictions (allowing only at most one particle on each site), 
deep in the species coexistence phase, i.e., for large predation rates, corresponding to a stable 
focal mean-field fixed point (stability matrix eigenvalues with negative real and non-vanishing 
imaginary parts). Both in the absence and presence of local density limitations, the coexistence 
phase is governed by remarkably strong spatio-temporal fluctuations: Striking spreading activity 
waves of prey closely followed by predators periodically sweep the system; any small surviving 
clusters of prey subsequently serve as sources for resurgent expanding prey-predator fronts [9j. 
An average over these weakly coupled local oscillations then yields the total population time 
traces depicted in figure [21 These spreading activity fronts appear especially sharp for the 
site-restricted model variants, as displayed in figure [U whereas in realizations with arbitrarily 
many particles per site, the fronts look more diffuse [23]. In either situation, one can employ 
stationary-state correlation functions to measure the spatial width ~ 10 ... 20 lattice sites of 
the spreading activity regions. At roughly the same length scale, the cross-correlations of the 
A and B particles peak at a positive value before slowly decaying to zero; at shorter distances, 
the prey are naturally anti-correlated with the predators [HI [23] . 




Figure 4. Snapshots illustrating the evolution in time (from left to right) of a two-dimensional stochastic 
lattice Lotka-Volterra model (with 512 x 512 sites), incorporating local occupation number restrictions in the 
species coexistence phase, with rates a — 4.0, /i — 0.1, A = 2.2, and D = 0, when the fixed point is a focus; the 
initial densities are a(0) = 1/3 — 6(0). The red, blue, and black dots respectively represent the prey, predators, 
and empty lattice sites. (Reproduced with permission from Ref. '[9,, p. 464.) 

In the absence of spatial degrees of freedom, the observed persistent population oscillations 
can be mathematically understood by performing a systematic van-Kampen expansion about 
the absorbing steady state |18j . The fluctuation corrections may then essentially be described 
through a damped harmonic oscillator driven by white noise that will on occasion resonantly 
incite large-amplitude excursions away from the stable fixed point in the phase plane. In 
our spatial systems, we may also interpret the persistent population oscillations in the species 
coexistence regime through a similar mechanism, as suggested by the bottom spiraling trajectory 
(for large predation rate A = 1.0) in the phase portrait depicted in figure [5l which was obtained 
in simulation runs for stochastic lattice Lotka-Volterra models with restricted site occupancy 
[9]. As the rate A is reduced (with all other parameters held constant) and the predators become 
less efficient, the stochastic lattice systems with site occupation restrictions qualitatively display 
the same scenarios as revealed by the mean- field analysis for eqs. dH) with finite prey carrying 
capacity: First, the focal stationary points in the phase plane are replaced by stable nodes 
(real stability matrix eigenvalues); the population oscillations then cease, and no interesting 
spatial structures aside from localized activity clusters with meek fluctuations are seen (c.f. the 
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Figure 5. Typical trajectories in the predator-prey phase space for a 512 x 512 stochastic lattice Lotka-Volterra 
system all initialized with a(0) = 1/3 = 6(0) and fixed rates a — 4.0, — 0.1, and D = 0, but different predation 
rates A = 0.15, 0.20, 0.40, 1.0. For small values of A (typically A < 0.4) the fixed point is a stable node, whereas 
for higher values of A one observe the characteristic spirals that indicate a focus in phase space. (Reproduced 
with permission from Ref. [9], p. 463.) 

trajectories for A = 0.4 and 0.2 in figure [5|). At a sufficiently small critical value Ac (~ 0.1688 
here, see figure [6]), a predator extinction threshold is encountered, and for A < Ac ultimately the 
prey population fills the entire lattice. 

For the predator population, the extinction threshold in stochastic spatial Lotka-Volterra 
models with local particle density restrictions represents a genuine continuous non-equilibrium 
phase transition in the thermodynamic and infinite-time limit. Since no conserved quantities or 
disorder are present, one expects this active to absorbing state phase transition to be described 
by the scaling exponents of critical directed percolation |27j-|32j. Heuristically, one may reason 
as follows: The prey density is essentially uniform and constant 6 ~ p near the critical point. 
The Lotka-Volterra reactions ([1]) then basically reduce to ^ — >• and A — t- A + A; but since 
the A population cannot multiply to arbitrarily large density (due to prey depletion), we need 
to add a growth-limiting reaction such as A + A ^ A, whereupon we arrive at the simplest 
microscopic reaction-diffusion model realization for the directed-percolation universality class 
[29t [3T| \32\ [37] . This assertion is indeed supported by careful analysis of Monte Carlo simulation 
data [Z]-IS], [II]-[I1], [inillZ]- We performed dynamical Monte Carlo simulations starting from 
a single active site with a predator particle in a lattice otherwise filled with prey, choosing 
reaction rates in the vicinity of the extinction threshold. The survival probability of predators 
at criticality is expected to decay algebraically as P{t) ~ , while the number of active sites 
with predators should grow according to the power law N{t) ~ [2HlEIl[32]i with 6' ^ 0.451 and 
9 ~ 0.230 for directed percolation in two dimensions |29[ I31|. Figure E] shows the corresponding 
effective exponents as functions of inverse time as measured in our Monte Carlo simulations 
for various values of A with a = 4.0, fj, = 0.1, and D = held fixed [9]. From these data we 
infer Ac ~ 0.1688 as best estimate for the critical predation rate (compare with the mean-field 
prediction Ac = = 0.1), and the extrapolation to t — )• oo yields very good agreement of the 
asymptotic critical exponents with the accepted directed-percolation values. 





Figure 6. Dynamical Monte Carlo simulation data to estimate the critical point and scaling exponents for the 
predator extinction threshold of the two-dimensional stochastic lattice Lotka-Volterra model with site occupation 
restrictions (on a 512 x 512 lattice): The effective scaling exponents —S'{t) vs. 1/t (left) and d{t) vs. 1/t (right) 
are depicted for four values of A (from top to bottom): 0.1690, 0.1689, 0.1688, and 0.1687, at fixed a = 4.0, 
H = 0.1, and D = 0. (Reproduced with permission from Ref. (£,, p. 470.) 



Simulations in one spatial dimension (on a circular domain) yield remarkable differences 
between model variants with and without site occupation number restrictions: In the former 
situation, the A and B particles quickly segregate into distinct domains, with the predation 
reactions occurring only at the boundary. The subsequent time evolution is governed by very 
slow coarsening induced by merging predator domains [9j. Without site occupation restrictions, 
in contrast we always observe the active coexistence state All the above statements of 

course pertain to sufficiently large lattices. In principle, any finite system with an absorbing 
steady state will eventually terminate in it; however, the associated survival times are expected 
to grow with system size according to a power law [19], and for our lattices are much longer 
than the duration of the simulations. 



3. Field-theoretic analysis 

3.1. Field theory representation 

In the remainder of this paper, I shall describe how stochastic fluctuations, internal reaction 
noise, and emerging correlations in spatial predator-prey models can be systematically captured 
by means of a field-theoretic representation of the associated classical master equation, which 
is then amenable to analytic approximations. Since the two-species Lotka-Volterra model ([1]) is 
defined via a diffusion-limited stochastic reaction system, we may employ the by now standard 
Doi-Peliti framework to map the associated master equation onto a field theory action |33]-|37j. 
This approach is based on the fact that at any time the configurations in such systems are 
completely enumerated through specifying the occupation numbers of each species per lattice 
site, and that all occurring stochastic processes merely modify these local integer occupation 
numbers. It is therefore natural to use bosonic creation and annihilation operators to formally 
represent the system's temporal evolution which is given in terms of a stochastic master equation. 
Subsequently the continuum limit can be taken, which in the time domain is most conveniently 
accomplished through a coherent-state path integral representation for the evolution operator. 

For the diffusion-limited reactions ([T]) in d spatial dimensions one thus arrives at the following 
action [S] (see also Ref. [3Tj) 



S[a,a;b,b] = i (fx i dt 



d(dt- DAV^)a + b (dt - Db 



+^{d — 1) a + (T il — b\ bb + Xib — d \ dab 



(5) 



with providing the statistical weight for any observables that must be functions of the 
local 'density' fields a{x,t) and b{x,t). The top line here obviously accounts for nearest- 
neighbor hopping processes in the continuum limit (through the inverse diffusion propagators 
with diffusivities and Db for the predators and prey, respectively). The bottom line in 
([5]) contains the stochastic reactions: spontaneous predator death with rate fi, prey birth with 
rate a, and predation with rate A. Note that each reaction process is represented by two 
contributions, originating from the gain and loss terms in the master equation. One can easily 
reconstruct these contributions in the action by noting that the second one directly reflects 
the reaction process itself through the annihilation operators a, b and creation operators d, 6, 
whereas the first one encodes the 'order' of the corresponding reaction (i.e., which powers of the 
concentrations a a and bb enter the rate equations). It is important to realize that the Doi-Peliti 
action faithfully contains all stochastic fluctuations associated with the underlying microscopic 
processes, namely discrete flnite- number fluctuations and internal reaction noise [57]. Following 
van Wijland's analysis [38J, restricted site occupation numbers or flnite local carrying capacities 
p for the prey species can be incorporated in this bosonic formalism through the replacement 
the B particle reproduction term. (Alternatively, a growth-limiting reaction 
such as B + B ^ B could have been added.) 

The associated classical fleld equations follow from the stationarity conditions 5S/5a = = 
5S/5b, always solved by d = 1 = 6 (actually just reflecting probability conservation ^37j), and 
6S/6a{x,t) = = 5S/6b{x,t), which yields precisely the mean-field rate equations augmented 
by diffusion terms. Indeed, for p = oo one arrives at eqs. ([2]), while expanding to first order 
in p'^^ recovers eq. dH. It is then convenient to perform a field shift according to d = 1 + d, 
b = 1 + 6, whereupon the action becomes, again to lowest order in the inverse carrying capacity. 



S[a, a; b,b] = J <fx j dt a(dt - DaV^ + p) a + b (dt -DeV^-a^b 

-ab^b + a p-^ {1 + bf bb^ - \{l + d){d -b)ab 



(6) 



In the following, the 'microscopic' field theory action Q will serve as the starting point (i) for 
further manipulations to identify the universality class of the continuous active to absorbing 
state phase transition at the predator extinction threshold, and (ii) to compute the fiuctuation- 
induced renormalization to lowest order in the predation rate for the population oscillation 
frequency and damping, as well as the diffusion coefficient in the two-species coexistence phase. 

3.2. Extinction transition and directed percolation 

Our goal is to construct an effective field theory ^ that describes the universal scaling properties 
near the non-equilibrium phase transition at Ac ~ p/p where the predators go extinct, and the 
prey fill the entire lattice: Og = 0, 6^ ~ p. Consequently we transform the action ([6]) to new 
fluctuating flelds c = bs — b with (c) = 0, and c = —b: 

S[d,a;c,c] = Jd'^xjdt d(dt - DaV'^ + p - \bs) a + c(dt - DbV'^ + {2bs/ p - I) c 

+a bs{2bs/p -l)c^ -a p'^ b^c^ -a {4hs/p - 1) c 
-a p~^ (1 + c^) c + 2fT p^^ c^ (c + bsc)c 



-A bs {d^ + (1 + d) a + A (1 + d) (d + c) a c 



(7) 



Next we note that the birth rate is a relevant parameter in the renormalization group sense, 
which scales to inflnity under scale transformations; this observation simply expresses the fact 



that fluctuations of the nearly uniform prey population become strongly suppressed through the 
'mass' term oc a for the c flelds. It is therefore appropriate to introduce rescaled fields (p = \fa c 
and (\) = \foc^ and subsequently take the limit a — )■ oo, which yields the drastically reduced 
effective action 



5oo [a, a; </>] = y d^x J dt d {dt — Da V"^ + fi — Xbg'j a — Xbg d'^ a + (f> (j) + bg (j)^ 



(8) 



As a final step, one needs to add a growth-limiting process for the predator population, for 
example through the binary coagulation reaction A + A ^ A with rate r. Since the fields (j) and 
(p only appear as a bilinear form in the action ([8]), they can readily be integrated out, leaving 



d'^x / dt 



i){dt + Da {ta - V^)^ V - uV (-0 - V') V' + V'^ 



(9) 



where ip = a^/Xb^pr, tp = d ^Jr/Xbg, rA = {fJ- — Xbs)/DA, and u = \/t Xbg. This new 
effective non-linear coupling u becomes dimensionless at dc = 4, signifying the upper critical 
dimension for this field theory. Near four dimensions, the quartic term oc r constitutes an 
irrelevant contribution in the renormalization group sense and may be omitted for the analysis 
of universal asymptotic power laws at the phase transition. The action ([9|) then becomes identical 
to Reggeon field theory, which is known to describe the critical scaling exponents for directed 
percolation [271 1321 EHl SO]- This mapping to Reggeon field theory [9] firmly corroborates 
the expectation that the predator extinction threshold is governed by the directed-percolation 
universality class [TIIH], [ISKTT], which features quite prominently in phase transitions 

to absorbing states |27u28] . even in multi-species systems [30j. The universal scaling properties of 
critical directed percolation are well-understood and quantitatively characterized to remarkable 
accuracy, both numerically through extensive Monte Carlo simulations and analytically by means 
of renormalization group calculations (for overviews, see Refs. [29 | I3H [32]). 



3.3. Fluctuation corrections in the coexistence phase 

In order to address fluctuation corrections in the predator-prey coexistence phase |42] . we start 
again from the Doi-Peliti field theory action ([6]), and introduce the proper fiuctuating fields 
c = a — {a) and d = b — (b): 



1 



pX 



+ AA+C 



(1 + B^) + d 



(10) 



Here, the mean-field values for the stationary densities have been taken into account already, 
such that the counter-terms A^ and B^^ which are naturally determined by the conditions 
(c) = = (d), contain only fluctuation contributions. The bilinear terms in the ensuing action 
may then readily be diagonalized by introducing new fields (p± and 
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with the mean-field (or 'bare') oscillation frequency and damping constant (see also Ref. j41j 



(11) 
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(12) 



Note that ujq = fia and 70 — >■ as p — )■ oo: There is no damping of the mean- field oscillations 
in the absence of local carrying capacity restrictions. 

In the following, we shall consider equal diffusivities Da = Dq = Db', the harmonic 
propagators in the diagonalized theory then read in Fourier space 



{4>±{q,uj)(t>±{q,uj'))o 



-iuj + Doq^ ± iuo + 70 



(13) 



Along with two two-point noise sources and several non-linear vertices, these propagators 
form the building blocks for the Feynman diagrams that graphically represent the different 
contributions in a perturbation expansion in terms of the non-linear coupling A [l2j. To lowest 
non-trivial ('one-loop') order, only the noise and three-point vertices are needed to determine 
the counter-terms Ac and Be, as well as to compute the fluctuation corrections to the bare 
propagators (jl3p . From the ensuing one-loop expressions, one may infer renormalized versions 
of the diffusivity Dr, oscillation frequency lor, and damping 7^. In addition, one finds that in 
the absence of site occupation restrictions (i.e., for infinite local prey carrying capacity p), the 
stochastic spatial fluctuations generate a damping term, just as seen in the lattice simulations. 
These perturbational calculations are fairly straightforward, but lengthy and somewhat tedious; 
details will be reported elsewhere [32]. Here I merely provide the explicit results for the 
renormalized parameters in several space dimensions. 

For d = 1 and d = 2, the expressions for the renormalized oscillation frequency become 
singular in the limit 70 — ?• 0; in the list below, only the leading terms in 70 are retained: 



d= 1 



d = 2: 



Dr 
IR = 

Dr 

IR = 




(14) 



(15) 



Notice that the infrared singularities encountered in the limit 70 — )• cancel for the renormalized 
diffusivity Dr and the fluctuation-generated damping ^r. In dimensions d < 2, the 



acquirmg 



leading fluctuation correction to the oscillation frequency diverges as (^0/70)^"'^''^, 
a logarithmic dependence in two dimensions; it is negative, and symmetric under formal rate 
exchange /i o a (c.f. the top lines in the above one-loop results for ujr). If we interpret 
7o in the above equations as a small, self-consistently determined damping, these features 
are in remarkable agreement with our earlier Monte Carlo observations displayed in figure [S] 
Fluctuations and correlations induced by the stochastic reaction processes induce a strong 
downward numerical renormalization of the oscillation frequency, with very similar functional 
dependence on the rates /i and a. 



In three dimensions, we may set the bare damping constant to zero (or p — > oo) to obtain 



d = 3 : Dr = D, 



(j) 4 V/i a 



LOR = 100 + 




128v/27r \DqJ 



^ 13 19 //I 13 /a iJL 

4 \l fj, A \ a 4: \fi a 



+ 0(A2), (16) 
hO(A2) . 



In higher dimensions d > 4, the fluctuation corrections become formally ultraviolet-divergent, 
and thus a finite cut-off A in momentum space must be implemented; e.g., in four dimensions 
one finds 
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- + - 
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Wij = Wo + 



256 TT \Do 



2 r 



. [1^. I. A4 

1-- J- In 1 + ^—2 

TT V CJ \ Wg/iJ^ 





+ 0(A") . 



(17) 



We finally remark that the effective expansion parameter in this fluctuation perturbation series 
in d dimensions is (A/wq) (ooo/Doy/^. 



4. Concluding remarks 

In conclusion, in this contribution I have reviewed the most striking features of stochastic 
predator-prey models on regular lattices that in the well-mixed mean-field limit reduce to 
the celebrated Lotka-Volterra model. It turns out that the spatially extended stochastic 
systems display both richer behavior than the associated deterministic rate equations, and 
are actually also more robust with respect to modifications of model and algorithmic details: 
Spatial predator-prey systems in the species coexistence phase are generically characterized 
by the emergence of persistent spatio-temporal structures, namely continually expanding and 
merging activity fronts, leading to transient oscillations for the total (or mean) particle densities. 
Fluctuations in the two-species coexistence phase are remarkably and unusually strong; they 
markedly alter the oscillation frequency as compared to the (linearized) mean-field prediction, 
and in addition generate damping. Restricting the (local) prey population through a growth- 
limiting finite carrying capacity induces a genuine continuous non-equilibrium extinction phase 
transition for the predators. I have also outlined how the Doi-Peliti field theory representation of 
the associated master equation can be employed to (i) demonstrate that this active to absorbing 
state transition is governed by the universal scaling exponents of critical directed percolation, 
and (ii) permits a systematic perturbational approach to compute the fluctuation-induced 
renormalizations of the population spreading and oscillation parameters in the coexistence phase. 

It remains to be elucidated which of the standard mathematical models in ecology, population 
dynamics, and chemical kinetics, many of which are frequently just discussed on the level 
of mean-field rate equations, are similarly strongly affected by stochastic fluctuations and 
intrinsic correlations. Perhaps unexpectedly, stochastic spatial variants of cyclic three-species 



predator-prey systems that are often referred to as rock-paper-scissors models represent an 
intriguing counter-example: Lattice simulations of these reaction-diffusion systems hardly show 
any noticeable fluctuation effects, both for model variants with conserved and non-conserved 
total particle number, despite the formation of striking spiral structures in the latter, so-called 
May-Leonard model, see Refs. [4:3\ 144] (and further references therein). 
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